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(-y-j ■ A simulated annealing algorithm is applied to the reconstruction 

of two-dimensional porous media with prescribed correlation func- 
tions. The experimental correlation function of an isotropic sample 
fj ' of Fontainebleau sandstone and a synthetic correlation function with 

C/3 \ damped oscillations are used in the reconstructions. To reduce the 

numerical effort we follow a proposal suggesting to evaluate the corre- 
lation functions only along certain directions. The results show, that 
this simplification yields significantly different micro-structures as com- 
pared to a full evaluation of the correlation function. In particular we 
find that the simplified reconstruction method introduces an artificial 
anisotropy that is originally not present. 
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I. INTRODUCTION 

A better understanding of the transport properties of random media, such as 
fluid flow in sandstones or electrical conductivity of composites requires the micro- 
structure as input ||l|-^. Digitized, three-dimensional micro-structures of natural 
sandstones are difficult and expensive to obtain j^. Thus there is a need for sim- 
ulation algorithms that are able to provide representative micro-structures from 
statistical probability functions. 

Recently various algorithms have been proposed for the reconstruction of ran- 
dom micro-structures [pHll[. In this paper we investigate a simulated annealing 
method that enforces agreement between the correlation functions of the original 
and the reconstructed micro-structure [pi. To save computation time the authors 
of g have evaluated the correlation functions only in certain directions assuming 
isotropy of the medium. The objective of this paper is to study the effect of this 
simplification on the final, reconstructed configurations. We test the effects on two 
examples: (i) the correlation function of a Fontainebleau sandstone and (ii) an ar- 
tificial correlation function with damped oscillations. The results show, that at 
least for the oscillating correlation function this yields significantly different config- 
urations. More importantly, the reconstructions are anisotropic as a result of the 
simplified evaluation of the correlation functions. 



II. THE RECONSTRUCTION METHOD 

We follow the reconstruction algorithm proposed in |8|]. The reconstruction is 
performed on a d-dimensional hyper-cubic lattice W/^ {d = 2 for the results presented 
below). Whether a lattice point lies within pore space or matrix space is indicated 
by the characteristic function 

(-\ I . fO for feP ._,. 

X{x) = x{^i,x2,....xa) = iy^ for f€M (^) 

with Xi = 0, 1, . . . , Mi — 1 where P denotes the pore space and M the matrix space 
of a two phase porous medium. The Xi are in units of the lattice spacing a. The 

porosity (}> is given as = jjYl,j=i (l ^ x{xj)] where N = 0^=1 -^^j i^ the total 
number of lattice sites. 

Simulated annealing is an iterative technique for combinatorial optimization prob- 
lems. The iteration steps are denoted by a subscript t. The optimum is found by 
lowering a fictitious temperature Tt that controls the acceptance or rejection of con- 
figurations with "energy" (or cost function) Et- The energy function used in our 
simulations is defined as 

Et^iZ ^^- E (•9'(^^) - aU^f (2) 

fc=l reDfc 

where Bj; C Z'' is a subset of lattice points and g^ is the fcth function of a set of 
J statistical probability functions calculated for the configuration of step t. For 
example g^ may be a /c-point correlation function. The real valued factor Wk is a 
weight for the fcth function. Hence, the energy can be understood as a measure for 
the deviations of the probability functions g^ from predefined reference functions 

5r^ef 

The simulated annealing algorithm consists of the following steps. 



1. Initialization: The O's and I's are randomly distributed with given porosity 

2. Two lattice points of different phase are chosen at random and exchanged. In 
this way the porosity (j) is conserved. 

3. The "energy" Et of the current configuration is calculated according to Equa- 
tion (1). 

4. The "temperature" Tt is adjusted according to a fixed cooling schedule. 

5. The new configuration created by the exchange of the two points is accepted 
with probability 

p = minh,expf ^ *~M j ■ (3) 

In case of rejection, the two points are restored and the old configuration is 
left unchanged. 

6. Return to step 2. 

As can be seen from Equation (H), configurations with lower energy are immedi- 
ately accepted while the acceptance of a configuration with higher energy is con- 
trolled by the temperature T. In order to obtain a configuration with minimal 
energy E oi, in other words, a configuration with minimal deviations of the proba- 
bility functions g'' from their reference functions g^^f, the temperature T has to be 
lowered in a suitable way. 

The algorithm is applicable to various functions g^. Most authors propose the 
two point correlation function |H,nO,nl| . Assuming homogeneity and ergodicity, the 
two-point correlation function can be defined for our case as 

,^ {X{x)x{x + r^) - {l - <P? 

9(-^ ^zr^2 • (4) 

where the average (. . .) indicates a spatial average over all lattice sites x. If gt-i{r) 
is known a single update step of the annealing process requires the recalculation of 
the correlations of the two pixels that are exchanged in step 2 with all other pixels. 
The numerical effort to obtain gt{r) is therefore proportional to N where A^ is the 
total number of lattice points. In the reconstruction of three-dimensional porous 
media with A^ ss 10^ . . . 10^ this leads to unacceptably long calculation times and 
therefore one has to find ways to speed up this calculation. 

One possibility for reducing the numerical effort is to truncate g at a value Tc for 
which g(r) pa Q for r — \r\ > r^ holds. Below we set g{r) = for |r| > r^ where 
Tc is a parameter in the reconstruction. For isotropic media, where g{r) — g(r) 
with \r\ = r, it was suggested in [pj to calculate g{r) only in certain directions 
by setting r ~ rck where el- is an arbitrary unit vector. In the two-dimensional 
reconstructions presented below Sk will be set to the radial unit vector in a polar 
coordinate system e/c = e^p^ = ex cos tpu + Cj, sin ^pk where Cx and Cj, are the unit 
vectors of the Cartesian coordinate system and ipk is the angle between Sx and e^p^. . 
Hence, instead of Equation (Eh we use 

fc, , {x{x)x{x + reu)) ' [l ~ (j)f 
^ (^) - J^^ (^) 

with r = 0, 1, . . . , Tc in the simplified reconstruction scheme. Since (^ is a set of J 
one-dimensional correlation functions the numerical effort is reduced by a factor of 
roughly Jrc/N as compared to (0). 



The above algorithm is now used to reconstruct two-dimensional, isotropic media 
with correlation function 

grc{{i~) = exp [~z) cos(tjr) (6) 

with r in units of the lattice spacing and u> = 1. The same function was used by 
the authors of ||] to exemplify their algorithm. In the evaluation of the correla- 
tion functions periodic boundary conditions are assumed. We use an exponential 
decrease of the temperature 

The remaining parameters are the lattice size Mi = M2 — 400 and the porosity 
(j) = 0.5. They are the same as in [p|. The following reconstructions have all 
been initialized with the same random seed. The algorithm terminated when the 
configuration did not change for 20000 subsequent update steps. 



III. RESULTS 

In the first reconstruction, the correlation function was calculated only in the 
horizontal and vertical direction, i.e. Equations (0) and dS) were used with J — 2, 
r\ = reo and 7^2 = ?'e,r/2- Both correlation functions g^^g^ had the same weight 
w\ = w-2, = 0.5 and the same reference function g\^^ = g^^f given in Equation (0) 
was used. The correlation function was truncated at Vc = 100. Figure 1 shows the 
final configuration of the reconstruction. A similar pattern was found in S . The 
pattern consists of stripes in direction of 6^/4 and in direction e_^/4. The distance 
between the stripes is determined by the cosine term. On a larger length scale, the 
pattern organizes into several regions in which all lines are parallel. The typical 
size of these regions is given by roughly 20 in view of Equation (|6|) and Figure y. 
Clearly, the pattern is not isotropic as it should be. The stripes are preferentially 
directed along the directions 6^/4 and e_^/4. Also, one expects that the oscillation 
frequency uj of the correlation function in direction e±7r/4 is not w = 1 but w = v2- 

Figure g shows the correlation functions for the configuration of Figure ^ in the 
directions of the unit vectors eo, 6^/2 and e^/4. The first and second have been used 
for the reconstruction and hence show good agreement with the reference function 
while the latter deviates drastically. The correlation in direction eV/4 is better 
described by the function 

fir) = i (exp (-0 cos (v^r) + exp (-0) . (8) 

shown as the dotted line. This may be interpreted as the arithmetic mean of the 
correlation function for regions (described above) with stripes perpendicular to the 
direction 6^/4 given by the first term in Equation (0) and the correlation function 
for regions with stripes parallel to the direction 67^/4 given by the second term. Of 
course, the same correlation function is found in direction e_T^ u. 

One step towards an isotropic reconstruction may be to use J = 4, and to force 
agreement of the correlations not only in two but in four directions eo, 6^/2, 6^/4 
and e_7r/4- The resulting configuration is shown in Figure |[ The pattern is signifi- 
cantly different from the pattern of Figure |l|. There are not only stripes parallel to 
the diagonal directions but more rounded formations and concentric circles at some 
points. The correlation functions for Figure are plotted in Figure 0. Surprisingly 



we find that it is not possible to obtain agreement with the reference correlation 
function. The resulting correlation functions in the horizontal, vertical and diago- 
nal direction all deviate strongly from the reference function especially at the first 
minimum. Additional simulations with different cooling schedules, damping factors 
and frequency of the correlation function did not give better agreement. We have 
also varied the system sizes from 200 x 200 up to 1000 x 1000. The results were 
identical. 

Figure [^ shows the two-dimensional reconstructions of a Fontainebleau sandstone 
with porosity = 0.135. The correlation function was truncated at re — 50. The 
reconstruction of Figure pk used the correlation function only in vertical eo and 
horizontal 6^/2 direction as suggested in |^]. The reconstruction shown in Fig- 
ure 5b, however, calculates the full two-dimensional correlation function according 
to ( 1 ) where the two-dimensional correlation function was radially binned in the 
calculations to obtain a one-dimensional function which can be compared to the 
one-dimensional, isotropic correlation function of the original sandstone. We em- 
phasize that in this calculation the two-dimensional correlation function was evalu- 
ated without restrictions or simplifications. Therefore, the calculation needed about 
50 times longer than the simplified reconstruction. Again, there are differences vis- 
ible although they are smaller than those in the reconstructions using (0). The 
shape of the pores in the full reconstruction (Figure gb) appear to be smoother 
and there is not as much "dust" visible as in the simplified (restricted) reconstruc- 
tion shown in Figure g|a. The correlation functions plotted in Figure ^ reveal also 
that the reconstructed micro-structure in Figure oa is strongly anisotropic while the 
micro-structure in Figure |5|b is isotropic as it should be. 

In summary we note that the statistical reconstruction of porous media with pre- 
defined two-point correlation function often requires some reduction of numerical 
effort. Especially the reconstruction of three-dimensional porous media in reason- 
able time does not seem possible without such simplifications in the calculation of 
the correlation function. Of course this problem is exacerbated if one wishes to 
include three-point or even higher order correlation functions. We applied a sim- 
plification proposed and used in [pj which samples the correlation function only in 
certain directions. The effect of this simplification on the final configurations may 
in some cases be negligible but in general configurations strong anisotropy and pat- 
terns which are significantly different from those of a proper isotropic reconstruction 
may appear as a result. 
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Figure 1: 
Figure 2: 



Figure 3: 
Figure 4: 
Figure 5: 



Figure 6: 



Simplified reconstruction of tlie damped oscillating correlation 
function given in Equation m) by restricting the correlation func- 
tion evaluation to the horizontal and vertical directions. 
Correlation functions for the reconstruction of Figure n^. The solid 
line is the reference function of the reconstruction given in Equation 
(||). + and X are the values of the correlation function of Figure ^ 
along the horizontal and vertical direction, respectively. The values 
of the correlation function in direction eV/4 areplotted with o. The 
dotted line is the estimate given in Equation ^ for the correlation 
function along the e^/4-direction. 

Simplified reconstruction of the damped oscillating correlation 
function of Equation (o) by restricting the correlation function eval- 



uation to the four eo, e^ 



72; ^TT/i 



and e_ 



-7r/4 



directions. 



Correlation functions for the reconstruction of Figure H. The ref- 
erence function is plotted as solid line. Note the mismatch in the 
first minimum and maximum. 

(a) Two-dimensional simplified reconstruction of the correlation 
function of a Fontainebleau sandstone by restricting the correlation 
function evaluation to the horizontal and vertical directions, (b) 
Reconstruction with the same reference correlation function as in 
(a) but complete evaluation of the correlation function, i.e. without 
directional restrictions. 

Correlation function for the reconstructions of Figure |^. The ref- 
erence function is plotted as solid line. The correlation functions 
calculated in horizontal, vertical and diagonal direction refer to the 
configuration of Figure pk. The complete correlation function is 
calculated from the configuration of Figure ||b. 
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Figure |l] 
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Figure H 
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Figure H 
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Figure kJ 
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Figure g 
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Figure H 
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